      SUBROUTINE LSHARM(RHO,TH,PH,F,NDIM,PNM,AJ,C,S)
COMPUTES LINE OF SIGHT FORCE FROM SPHERICAL HARMONICS.
      DOUBLE PRECISION PNM(NDIM,NDIM),PPNM,F(3),STH,CTH,CO,R1
     1,CPH,SPH,CT2,C1,C2,C3,C4,C5,S2,S1,EN,RRN,EM,S3,S4,S5,S3M,S4M,S5M,
     2PM,SMP,CMP,RR,AJ(NDIM),C(NDIM,NDIM),S(NDIM,NDIM),C6,C7,C8
      DOUBLE PRECISION COM(1)
      DOUBLE PRECISION TH,PH,RC,RHO,U
      COMMON /PRAM/ NMAS, NPTS, NDISK
      COMMON RC,U,COM
C     DOUBLE PRECISION POUT(16)
      DATA CO /.1745329251994430D-1/
C     DEFINE AJ(I)=COM(I+IA)
C     DEFINE C(I,K)=COM(K*NPTS+I+IA)
C     DEFINE S(I,K)=COM(NPTS*NPTS+K*NPTS+I+IA)
      IA = 4*NMAS
      RR=RC/RHO
      R1=1.D0/(RHO*RHO)
      CTH=DCOS(TH*CO)
      STH=DSIN(TH*CO)
      CPH=DCOS(PH*CO)
      SPH=DSIN(PH*CO)
      CT2=CTH*CTH
      C1=STH*CTH*CPH*U/RHO**2
      C2=CTH*CPH*U/RHO**2
      C3=-C1
      C4=SPH*U/(RHO*RHO*CTH)
      C5=-C2
      C6=STH*CTH*SPH*(U/RHO**2)
      C7=CTH*SPH*(U/RHO**2)
      C8=(-CPH/CTH)*(U/RHO**2)
      S2=0.D0
      S1=0.D0
      PNM(1,1)=1.D0
      PNM(2,1)=STH
      NP=NPTS+1
      DO 1 N=3,NP
      EN=DBLE(FLOAT(N-1))
      PNM(N,1)=STH*(2.D0*EN-1.D0)*PNM(N-1,1)/EN-PNM(N-2,1)*(EN-1.D0)/EN
      RRN=RR**(N-1)
      PPNM=(EN*PNM(N-1,1)-EN*STH*PNM(N,1))/CT2
      S1=S1+AJ(N-1)*PPNM*RRN
    1 S2=S2+(EN+1.D0)*AJ(N-1)*PNM(N,1)*RRN
      DO 2 M=2,NP
      EM=DBLE(FLOAT(M-1))
      L=M+2
      PNM(M,M)=(2.D0*EM-1.D0)*CTH*PNM(M-1,M-1)
      IF(M.EQ.NP)GO TO 2
      PNM(M+1,M)=(2.D0*EM+1.D0)*STH*PNM(M,M)
      IF(L.GT.NP)GO TO 2
      DO 3 N=L,NP
      EN=DBLE(FLOAT(N-1))
      PNM(N,M)=((2.D0*EN-1.D0)*STH*PNM(N-1,M)-(EN-1.D0+EM)*PNM(N-2,M))/
     1(EN-EM)
    3 CONTINUE
    2 CONTINUE
      S3=0.D0
      S4=0.D0
      S5=0.D0
      DO 4 N=2,NPTS
      S3M=0.D0
      S4M=0.D0
      S5M=0.D0
      NN=N+1
      RRN=RR**N
      EN=DBLE(FLOAT(N))
      DO 5 M=1,N
      MM=M+1
      EM=DBLE(FLOAT(M))
      PM=DMOD(EM*PH,360.D0)
      SMP=DSIN(PM*CO)
      CMP=DCOS(PM*CO)
      PPNM=(-EN*STH*PNM(NN,MM)+PNM(NN-1,MM)*(EN+EM))/CT2
      S3M=S3M+PPNM*(C(N,M)*CMP+S(N,M)*SMP)*RRN
      S4M=S4M+EM*PNM(NN,MM)*(C(N,M)*SMP-S(N,M)*CMP)*RRN
      S5M=S5M+(EN+1.D0)*PNM(NN,MM)*(C(N,M)*CMP+SMP*S(N,M))*RRN
    5 CONTINUE
      S3=S3+S3M
      S4=S4+S4M
      S5=S5+S5M
    4 CONTINUE
C     WRITE(6,6) S1,S2,S3,S4,S5
C     WRITE(6,7) C1,C2,C3,C4,C5,C6,C7,C8
C     WRITE(6,8) SPH,CPH,CTH,STH,U,RHO,TH,PH
C   6 FORMAT(1H ,'S=',5D16.8)
C   7 FORMAT(1H ,'C=',8D15.8)
C   8 FORMAT(1H ,'SPH,CPH,CTH,STH,U,RHO,THETA,PHI=',8D11.5)
      F(1)=C1*S1+C2*S2+C3*S3+C4*S4+C5*S5
      F(2)=C6*(S1-S3)+C7*(S2-S5)+C8*S4
      F(3)=U*(CTH*CTH*(-S1+S3)+STH*(S2-S5))/(RHO*RHO)
C     WRITE(6,9) F(1),F(2),F(3)
C   9 FORMAT(1H ,'F(1,2,3)=',3D16.8)
C     DO 11 IOUT=1,NPTS
C     DO 10 IOUT2=1,NPTS
C     POUT(IOUT2)=PNM(IOUT,IOUT2)
C  10 CONTINUE
CC    WRITE(6,12) IOUT,(POUT(MM),MM=1,NPTS)
C  11 CONTINUE
C 12  FORMAT(1H ,' PNM(',I3,',1)=',4D16.8)
      RETURN
      END
